Anchoring effects at the isotropic-nematic interface in liquid crystals 
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The isotropic-to-nematic transition in liquid crystals is studied in d = 3 spatial dimensions. A 
simulation method is proposed to measure the angle dependent interfacial tension y(8), with 9 the 
anchoring angle of the nematic phase at the interface. In addition, an alternative liquid crystal 
model is introduced, defined on a lattice. The advantage of the lattice model is that accurate 
simulations of anchoring effects become possible. For the lattice model, y(6) depends sensitively 
on the nearest-neighbor pair interaction, and both stable and metastable anchoring angles can be 
detected. We also measure "/(O) for an off-lattice fluid of soft rods. For soft rods, only one stable 
anchoring angle is found, corresponding to homogeneous alignment of the nematic director in the 
plane of the interface. This finding is in agreement with most theoretical predictions obtained for 
hard rods. 

PACS numbers: 83.80.Xz, 68.05.-n, 68.03.Cd, 64.70.Md, 61.30.Hn 



I. INTRODUCTION 

A fluid consisting of elongated molecules is more diffi- 
cult to describe than one in which the molecules are sim- 
ply spherical. In the case of elongated molecules, there 
are not only translational degrees of freedom, but also 
orientational ones. This additional complexity gives rise 
to many interesting effects, not found in spheres. For 
example, infinitely slender rods in three dimensions un- 
dergo a first-order phase transition from an isotropic to 
a nematic phase, provided the density is sufficiently high 
[l[ . Both the isotropic and the nematic phase lack trans- 
lational order, but in the nematic phase the rods have 
aligned, giving rise to long-range orientational order. 

The orientation of the nematic phase is an important 
quantity. In applications involving nematics at walls, the 
angle of the nematic director at the wall is often crucial. 
This angle is called the tilt or anchoring angle. Typi- 
cally, there is a preferred tilt angle the nematic phase 
will assume, but the precise value depends sensitively 
on factors such as surface chemistry, particle shape, and 
temperature @, [H, 0|. Similarly, anchoring effects also 
occur at the isotropic-to-nematic (IN) transition. The 
first-order nature of that transition implies phase coex- 
istence, whereby isotropic domains coexist with nematic 
domains, separated by interfaces. As Fig. [T] shows, the 




FIG. 1: Schematic representation of isotropic-nematic phase 
coexistence. The isotropic phase is on the left, the nematic 
on the right. Shown are (a) homogeneous anchoring, and (b) 
homeotropic anchoring. 



orientation of the nematic phase with respect to the in- 
terface becomes an additional parameter. In Fig. [lja), 
the nematic director points in the plane of the interface, 
which is called planar or homogeneous alignment. In 
Fig. [IJb), the director is perpendicular to the interface, 
which is known as homeotropic alignment. 

From symmetry considerations alone, it is clear that 
homogeneous and homeotropic anchoring are different. 
For homeotropic anchoring, there is still rotational sym- 
metry around the interface normal; for homogeneous an- 
choring, no such symmetry is present. This difference is 
known to affect the spectrum of capillary waves. For ho- 
mogeneous anchoring, the spectrum becomes anisotropic 
in the short wavelength limit [f| @, 0>§]. In contrast, for 
homeotropic anchoring, the spectrum remains isotropic 
at all wavelengths. In other words, as this example shows, 
the anchoring angle affects the interfacial properties qual- 
itatively. Given a set of particle interactions, it is there- 
fore important to be able to predict the anchoring angle. 
This has lead to the concept of an angle dependent in- 
terfacial tension 7(0), with 9 the tilt or anchoring angle. 
Here, 9 is defined as the angle between the nematic di- 
rector and the plane of the IN interface. Homogeneous 
(9 = 0) and homeotropic (9 = 90) anchoring are most 
common, although 9 could, in principle, be anywhere be- 
tween and 90 degrees. In theoretical investigations, the 
anchoring angle is given by the angle which minimizes 
j(9). For hard rods, this is typically 9 = 0, correspond- 
ing to homogeneous anchoring 

@, Ha HI HI mp, but 

the precise behavior is quite subtle. For example, the 
results of Ref. H also suggest that homeotropic anchoring 
may be metastable. In addition, for very short rods, an- 
choring angles between and 90 degrees have also been 
reported fl2| . 

Unfortunately, it remains difficult to verify these the- 
oretical findings in a computer simulation. On the one 
hand, efficient simulation methodology for problems of 
this kind is scarce. The state-of-the-art is to extract j(0 ) 
from the anisotropy of the pressure tensor 0, [lH, fig . 
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a technique which is somewhat prone to statistical er- 
ror. On the other hand, the particle interactions used 
in many theoretical investigations are not convenient for 
simulations. The hard rod potential, for instance, often 
used in theory, gives rise to a very small interfacial ten- 
sion. In order to stabilize the IN interface, simulations of 
hard rods require huge system sizes, implying long equi- 
libration times and, consequently, data with considerable 
statistical uncertainty. 

The purpose of this paper is to improve on this state 
of affairs. The primary aim is to present a simulation 
method capable of measuring the angle dependent in- 
terfacial tension accurately. The method is presented in 
Section|TTl As it turns out, the method is general, and ap- 
plies to lyotropic (density driven) systems, such as rods or 
platelets, as well as to thermotropic (temperature driven) 
lattice systems. The second aim is to introduce a new 
liquid crystal model, one which is easy to simulate, but 
which nevertheless features an IN transition with anchor- 
ing effects. The model we propose is defined on a lattice, 
and resembles the Lebwohl-Lasher (LL) model [13], but 
with two essential modifications. Since the model is easy 
to simulate, it lends itself perfectly for an investigation 
of anchoring effects. The liquid crystal model, and the 
subsequent determination of its 7(0), are presented in 
Section Um Next, in Section UVl we determine j(8) for a 
fluid of soft rods. These particles are already more com- 
plicated to simulate. Nevertheless, guided by the experi- 
ence obtained for the simple lattice model, a meaningful 
interpretation of the simulation data is possible. We end 
with a summary and outlook in the last section. 



II. SIMULATION METHOD 

In this Section, we present our method to extract the 
angle dependent interfacial tension "f(6) in liquid crys- 
tals. The use of so-called order parameters is crucial for 
our method: suitable order parameters are therefore dis- 
cussed first. Next, we show how the order parameter dis- 
tribution may be used to obtain phase coexistence prop- 
erties, as well as interfacial tensions, which summarizes 
the key ingredients of previous work [l8|, [h| . Finally, we 
show how this methodology can be modified, to also cap- 
ture the angular dependence of the interfacial tension, by 
means of a simple constraint. 



A. Order parameters 

Since we are dealing with the IN transition in liquid 
crystals, a suitable order parameter is the nematic order 
parameter S, defined as the maximum eigenvalue of the 
orientational tensor Q, whose elements read as: 

1 N 

Qap = 77^7 / J (Sdjgdjp - Sap) . (1) 
i=i 



Here, di a is the a component (a = x,y,z) of the orien- 
tation di of molecule i (normalized to unity) , 8 a p is the 
Kronecker delta, and N the number of molecules. Note 
that S is invariant under the inversion di — > — di of single 
molecules, which is the characteristic symmetry of liquid 
crystals, and also that S does not depend on the center 
of mass coordinates. In the isotropic phase, S is close 
to zero. In the nematic phase, where the molecules have 
aligned, S is close to unity. Another important quantity 
is the (normalized) eigenvector ft = (n x ,n y ,n z ) associ- 
ated with S. The vector n is called the director, and 
it corresponds to the overall preferred direction of the 
molecular orientations in the nematic phase. Again, the 
directions n — > — n are equivalent: the convention in this 
work is to pick the vector with n z > 0. 

The nematic order parameter, being zero in the 
isotropic phase and (close to) unity in the nematic phase, 
is a convenient quantity to detect the IN transition in 
liquid crystals. However, different quantities may be 
used as well. For example, in thermotropic (temperature 
driven) liquid crystals, there is also an energy difference 
between the isotropic (high energy) and nematic (low en- 
ergy) phase. Therefore, in thermotropic systems, energy 
may also be used as order parameter. Similarly, in ly- 
otropic (density driven) liquid crystals, such as studied 
by Onsager [l[ , there is also a density difference between 
the isotropic (low density) and nematic (high density) 
phase. Therefore, in lyotropic systems, density is also a 
valid order parameter. 



B. Order parameter distributions 

Our method to obtain 7(0) is based on the order pa- 
rameter distribution P(X), defined as the probability to 
observe the order parameter X during the simulation. 
For liquid crystals, suitable choices for X were given 
above. As is well known, at a first-order phase transition, 
the distribution P(X) becomes double-peaked (bimodal). 
An example is provided in Fig.J^a), which shows the en- 
ergy distribution P{E) of a thermotropic liquid crystal 
(details are provided in Sect ion [Hi]). In thermotropic sys- 
tems, P(E) becomes bimodal at the transition tempera- 
ture. The precise value is determined using the "equal- 
area" rule [20] . whereby the temperature is tuned such 
that the area under both peaks is equal. Of course, in 
a lyotropic system, one would need to tune the chemical 
potential. 

From the bimodal energy distribution of Fig. [^a), bulk 
properties can readily be extracted. The peak at high 
energy, for example, yields the energy density of the 
isotropic phase; the peak at low energy of the nematic 
phase. Even more information is contained in the loga- 
rithm W = lnP(A), see Fig. [^b). Note that W corre- 
sponds to minus the free energy of the system. We now 
observe a distinct flat region between the peaks. The 
origin of this flat region can be understood from simula- 
tion snapshots, shown schematically in Fig. [3] When the 
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FIG. 2: (a) Coexistence distribution P(E) of a thermotropic 
liquid crystal, interacting via Eq.@ with p = 10 and = 
0.5, at the transition (inverse) temperature e* ~ 1.188. The 
simulation box dimensions were L — 15 and D = 40. (b) The 
logarithm of the same distribution. 



in the free energy at all, and hence a fiat region in W. 

The presence of a flat region in W naturally allows for 
an estimate of the interfacial tension [2l[ ■ In these cases, 
the height of the free energy barrier AF in Fig. [2b) may 
be associated with the free energy cost of having two in- 
terfaces in the system. Since the interfacial tension is 
defined as the excess free energy per unit area, one sim- 
ply has 7 = AF/2A, with A the area of one interface. It 
was later recognized that [22| , in an elongated L x L x D 
simulation box, with D 3> L, the interfaces form perpen- 
dicular to the elongated direction, since this minimizes 
the total amount of interface in the system. This leads 
to A = L 2 , and consequently 

7 = AF/2L 2 . (2) 

In previous work, the above ideas were successfully ap- 
plied to the IN transition in fluids of rods [H, [l9[ and 
platelets (23[. Implementation details are also provided 
in these references. Of particular importance is the use 
of a biased sampling scheme (24j, such that the simula- 
tion frequently traverses between the isotropic and the 
nematic phase. 




FIG. 3: Schematic simulation snapshots in (a) the bulk 
isotropic phase, (b) the coexistence region, and (c) the bulk 
nematic phase. 



system is in the high-energy peak, simulation snapshots 
reveal a homogeneous isotropic phase (a). In the low- 
energy peak, snapshots reveal a homogeneous nematic 
phase (c). At intermediate energy, coexistence between 
an isotropic and nematic domain is revealed, separated 
by an interface (b). Note that, due to periodic bound- 
ary conditions, two such interfaces are actually present. 
Provided the simulation box is large enough so as to ac- 
commodate two non-interacting interfaces, the order pa- 
rameter can be varied (over a limited range) with no cost 



C. Measuring 7(6*) 

Next, we describe how to modify the above methodol- 
ogy to also extract the angular dependence of the inter- 
facial tension. We again use an elongated simulation box 
with periodic boundary conditions. The box is spanned 
by the vectors Lx, Ly, and Dz, with D » L. As usual, 
x = (1,0,0), y — (0,1,0), and z — (0,0,1) denote stan- 
dard Cartesian unit vectors. The key additional ingredi- 
ent is to add a constraint to the Hamiltonian, such that 
the total energy of the system becomes E — Eq + E c . 
Here, Eq is the energy of the unconstrained system. For 
example, in a thermotropic system, Eq could be the LL 
potential. In a lyotropic system, it could be the potential 
of hard rods. The constraint energy E c should fulfill two 
criteria: 

1. In the bulk isotropic and nematic phase, the influ- 
ence of the constraint must vanish. In other words, 
E c may not affect the bulk properties of the uncon- 
strained system. 

2. In the coexistence region, where the system 
schematically resembles Fig. G2b) , the director n of 
the nematic phase must point along some specified 
tilt angle 9. 

As it turns out, a suitable constraint can be written as: 

JO |90 — arccos \ri ■ z\ — 9 t \ < 5, 
^c= < (3) 
00 otherwise. 

Here, n is the nematic director, defined in Section Til Al 
The angles 9t and 8 are inputs of the method, and 
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must be specified beforehand. By using the constraint, 
only states whose angle between director and xy-j>\anc 
is within 9 t ± S are retained, while all other states are 
rejected. 

For large systems, Eq.([3]) does not affect bulk proper- 
ties, since bulk properties are insensitive to the overall 
orientation of the phase. In contrast, in the coexistence 
region, the constraint has a dramatic effect. In these 
cases, approximately half of the simulation box is filled 
with an isotropic domain, and the other half with a nc- 
matic domain, see Fig. [3Jb) . Due to the constraint, the 
angle between the director of the nematic domain and 
the xy-plane is within 9t ± S. At the same time, the use 
of an elongated simulation box forces the interfaces to 
form in the xy-plane as well. In other words, by setting 
9 t , the anchoring angle 9 can be fixed. More precisely, 
one has 9 — 9 t . Naturally, the threshold angle 5 should 
be chosen as small as possible, while, at the same time, 
maintaining reasonable simulational efficiency. The op- 
timal value is model dependent, and best obtained using 
trial-and-error. 

The idea to obtain 7(6*) is now clear. We first choose 
a tilt angle 9 of interest. Next, we measure the order 
parameter distribution P(X), in an elongated simulation 
box, using the methodology of Section III Bl In addition, 
we incorporate the constraint of Eq.([3]) in the simula- 
tions, using 8t = 9. The peak positions in P(X) should 
again yield the bulk properties of the coexisting isotropic 
and nematic phase. The barrier AF, see Fig. Hfb), can 
be plugged into Eq.([2]) to obtain the interfacial tension 
at the chosen anchoring angle 9. Since bulk properties 
should not be affected by the constraint, we expect the 
peak positions in P(X) to coincide with those of an un- 
constrained simulation. In contrast to the unconstrained 
simulations, a dependence of the interfacial tension on 
the anchoring angle 9 is anticipated. To what extent 
these expectations are met in actual simulations will be 
investigated next. 



III. RESULTS: LATTICE SIMULATIONS 

A. Lattice model and motivation 

As announced in the Introduction, we first test our 
method in a lattice model of a thermotropic liquid crys- 
tal. The aim is to measure "f(9). The simulations are 
performed on a three-dimensional periodic lattice of size 
V = L x L x D, with D ^> L. To each lattice site i, a 
liquid crystal is attached with (normalized) orientation 
di. The liquid crystals interact via the potential 



<i,j> 



di\ p , 



(4) 



(a) side-side alignment 



(b) head-head alignment 

FIG. 4: Illustration of the spatial anisotropy in the liquid 
crystal pair interaction. In a realistic system, the energies 
of the above two configurations will generally differ. The LL 
model, however, makes no distinction. 



the temperature, and fee the Boltzmann constant. The 
"anisotropy" parameter is given by 



= l + v (di - f i:j ) 2 + (dj ■ fijf 



(5) 



where the summation is over nearest neighbors, coupling 
constant e, and exponent p > 0. In what follows, factors 
of k^T are absorbed in the coupling constant e, with T 



with Tij a unit vector pointing from site i to j, and v a 
parameter between —0.5 < v < 0.5. On a cubical lattice, 
each site has six nearest neighbors. Consequently, there 
are only three possible axes along which the vectors 
can be oriented. 

For p = 2 and v = 0, this model reduces exactly to 
the LL model [l7|. In this case, the model exhibits a 
first-order IN transition, but it is very weak 25]. This 
makes the LL model rather inconvenient for our pur- 
poses. For example, to stabilize two interfaces so as to 
recover the coexistence of Fig. [3Jb), huge systems would 
be required. Such large-scale simulations are not the aim 
of the present work, and so we have chosen to modify 
the interactions appropriately. More precisely, we use 
a larger exponent in Eq.(U]), namely p = 10. The ef- 
fect of this is a sharper pair interaction, meaning that 
neighboring molecules only lower their energy when they 
are closely aligned. It is known that, under such inter- 
actions, first-order phase transitions become enhanced 
[H HE H, [U, HI (even in two dimensions). 

By using p = 10 in Eq.Q, the model is expected to ex- 
hibit a strong first-order IN transition. Nevertheless, this 
is not sufficient to study anchoring because, for v = 0, 
the interactions are spatially isotropic. In other words, 
the interactions do not depend on the relative positions 
of the molecules, and so they cannot produce any an- 
choring effects at the IN interface [3lJ. In realistic sys- 
tems, the particle interactions are typically anisotropic, 
see Fig. HJ Shown are two liquid crystal arrangements, 
labeled (a) and (b). Even though the orientations of the 
molecules are identical in both cases, it is clear that the 
energies need not be the same. In the LL model, how- 
ever, for which v = 0, there is no distinction between 
the two arrangements. In order to nevertheless study an- 
choring effects, we allow v ^ in Eq.([5]), in which case 
the model does make the distinction. More precisely, we 
have (Tij — 1 for case (a), and cr^ = 1 + 2^ for case (b). 
By choosing v < 0, side-side alignment is energetically 
favored; choosing v > favors head- head alignment. For 
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FIG. 5: Bulk properties of Eq.Q with p = 10 as a function of 
v. Points are actual simulation data; curves serve to guide the 
eye. (a) Variation of e* with v. (b) Binodal curves, showing 
the energy density E/V of the isotropic phase (top curve), 
and of the nematic phase (lower curve), as a function of v. 
At energy densities between the curves, coexistence between 
isotropic and nematic domains occurs; simulation snapshots 
will then schematically resemble Fig. [3jb). 



FIG. 6: Order parameter distributions of Eq.@ with p — 10 
and v — 0, obtained using box dimensions L = 10 and 
D — 30. Shown are the logarithm of the energy distribution 
(a), and of the nematic order parameter (b). In each case, 
three distributions are shown, corresponding to homogeneous 
and homeotropic enforced anchoring, as well as no enforced 
anchoring direction. The curves overlap almost perfectly, in- 
dicating the absence of any anchoring effects. 



v = 0, the interactions are isotropic, in which case no 
particular alignment is preferred. 



B. Bulk phase behavior 

We first determine the bulk behavior of Eq. (J3J) . Recall 
that we keep the exponent fixed at p = 10. The aim is to 
measure the variation of the bulk properties as a function 
of v. More precisely, we consider the transition inverse 
temperature e*, and the coexistence energy densities of 
the isotropic and nematic phase. To this end, we use 
the simulation methodology of Section III Bl without the 
constraint of Eq.©. The energy distribution P(E) is 
measured in a MC simulation, using a biased sampling 
scheme [24J, at e = 0. Histogram reweighting 32] is 
used to determine the value of e for which the "equal- 
area" rule is obeyed, yielding e*. The energy densities 
are then read-off from the peak positions. An example 
distribution P(E) is shown in Fig. [5] The simulations 
are performed using single particle MC moves, whereby 
a random orientation is assigned to a randomly selected 
lattice site, accepted with the Metropolis criterion (33|. 
Typical lattice sizes are L = 10 - 20 and D = 20 - 40. 
The CPU time required to obtain P{E) accurately for a 
large system is around 48 hours. 

The variation of e* with v is shown in Fig. [UJa). The 
behavior is simply monotonic: by increasing v, e* goes 
down. The energy densities, shown in Fig. EJb), reveal 
more interesting behavior. By decreasing the energy 



difference between the isotropic and the nematic phase 
becomes smaller. In other words, the transition becomes 
weaker. The simulations do not rule out that the curves 
meet when v becomes sufficiently negative, possibly ter- 
minating in a critical point, but clearly additional efforts 
are required to resolve this. All that matters for the 
present work, however, is the fact that Fig. [HJb) reveals 
a large coexistence region, over a substantial range of 
v values. This confirms our expectation that, by using 
p = 10 in Eq.(|l]), the first-order nature of the transition 
is enhanced significantly. This makes the model ideal to 
study anchoring effects, with which we proceed next. 



C. Anchoring effects for v — 

As a benchmark, we consider Eq.Q with v = 0, us- 
ing the newly proposed method. Note that, for v = 0, 
the model is spatially isotropic, and so we do not expect 
to see any anchoring effects. We MC simulate Eq.Q as 
before, with the constraint of Eq.Q explicitly included. 
Two anchoring conditions are considered: homogeneous 
and homeotropic. Recall that the anchoring is set via 6 t 
in Eq.©. For the threshold angle, we use 5 = 0.75 de- 
grees. In addition to P(E), we also measure P(S), with 
S the nematic order parameter defined in Section III Al 
For completeness, we mention that our simulations are 
performed using a bias on the nematic order parameter 
S, see details in Ref. 

The resulting energy distributions are given in 
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FIG. 7: Transition inverse temperature e* versus anchoring 
angle 9 for Eq.(j4| with p — 10 and three values of v as indi- 
cated. Closed symbols are actual simulation data; lines serve 
to guide the eye. The data were obtained using box dimen- 
sions L — 15 and D — 40. Note the much finer scale in (a) 
compared to (b) and (c). 



Fig. EJa), which actually shows three distributions. 
Shown are the two distributions obtained using the new 
method, corresponding to homogeneous and homeotropic 
anchoring, as well as the distribution obtained without 
any enforced anchoring. The striking feature is that the 
curves overlap almost perfectly. This result is crucial 
because it demonstrates the consistency of the method. 
First of all, the insensitivity of the peak positions with re- 
spect to the enforced anchoring, confirms that bulk prop- 
erties are not affected by the constraint. In addition, 
we find that the barrier AF, defined in Fig. [2fb), also 
does not depend on the anchoring condition. In other 
words, the interfacial tension is independent of the tilt 
angle, which is precisely what one expects for an isotropic 
potential. Additional confirmation of the consistency of 
the new method is provided in Fig. HJb), which shows 
the corresponding distributions In P(S) of the nematic 
order parameter. Again, the curves overlap almost per- 
fectly. Note also that, for the interfacial tension, it does 
not matter whether one reads-off the barrier height in 
In F(F) or In P(S). As Fig. [S] shows, the barriers are 
nearly equal (the slight variation gives an indication of 
the statistical uncertainty). 

We have repeated the above analysis using larger lat- 
tices, considering also tilt angles between and 90 de- 
grees. Shown in Fig. [TJa) is the transition inverse tem- 
perature e* versus 9. As expected, for the spatially 
isotropic case, e* is insensitive to 9, and we obtain 
e* = 1.5860 ± 0.0005. Shown in Fig. lis the angle de- 
pendent interfacial tension 7(0), as extracted from the 
barrier AF in \nP(S) and using Eq.@. Here, AF was 
taken to be the average height of the peaks, measured 
with respect to the flat region. As expected, the inter- 




FIG. 8: Angle dependent interfacial tension y(9) of Eq.Q us- 
ing p = 10 and three values of v as indicated. Closed symbols 
are actual simulation data; lines serve to guide the eye. The 
data were obtained using box dimensions L — 15 and D — 40. 
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FIG. 9: Logarithm of the nematic order parameter distribu- 
tion P(S), at coexistence, of Eq.Q with p = 10 and v — 0.5. 
Shown is In P(S) for various imposed anchoring angles 9, with 
9 the angle between the nematic director and the plane of the 
IN interface. The distributions were obtained using box di- 
mensions L = 15 and D — 40. 



facial tension does not display any pronounced 9 depen- 
dence (the variation stays below 2%). For v — 0, we thus 
find 7 = 0.108 ± 0.002 k-^T per squared lattice spacing, 
independent of the tilt angle. 



D. Anchoring effects for v = 0.5 

Having verified that the spatially isotropic case v = 
does not reveal any anchoring effects, we now consider 
Eq.flUl using v = 0.5. In this case, the model becomes 
anisotropic, and the tilt angle of the nematic phase with 
respect to the IN interface should become a relevant pa- 
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FIG. 10: Profiles S(z) (dashed curves) and a(z) (solid curves) of Eq.Q with p = 10 and v — 0.5, for various imposed tilt 
angles 6. The profiles were obtained at the transition inverse temperature e*, overall nematic order parameter S = 0.4, and 
box dimensions L = 40 and D = 100. 



rameter. For a number of tilt angles, we have measured 
the order parameter distribution In P(S) at coexistence; 
recall that the tilt angle is set via the constraint of Eq.Q. 
Typical distributions are plotted in Fig. [5] Shown in (a) 
are distributions for tilt angles close to homogeneous and 
homeotropic anchoring; shown in (b) are distributions for 
two "in-between" tilt angles. Also shown in (a) is the 
"unconstrained" distribution, which one obtains without 
imposing the constraint of Eq. ^ . For tilt angles that are 
close to 8 = or 8 = 90, the distributions behave as ex- 
pected: they are bimodal, and also exhibit a pronounced 
flat region between the peaks. In addition, we observe 
that the barrier height, defined in Fig. [3J depends pro- 
foundly on the imposed tilt angle. Since the barrier is re- 
lated to the interfacial tension, via Eq. (J2J) , we can already 
see that anchoring effects are present. Interestingly, for 
the "in-between" tilt angles, bimodal distributions can 
also be identified, but the region between the peaks is 
not quite flat, see the arrow in Fig. [5fb). This suggests 
that, for these "in-between" angles, the constraint does 
not quite produce the IN coexistence scenario of Fig. [3l 
but rather something else. 

To verify what is going on, we have generated a number 
of snapshots, at fixed nematic order parameter S = 0.4. 
For all distributions in Fig. [9l this value is well between 
the peak positions. The snapshots are generated at the 
transition inverse temperature e* using MC simulation. 
After equilibration, we collect the profiles S(z) and a(z). 
Here, S(z) is the nematic order parameter in the z-th 
LxL slab perpendicular to the elongated z-direction, and 
a(z) the angle between the director in that slab and the 
xy-plane. The profiles are shown in Fig. 1101 for the same 
tilt angles as studied in Fig. [9] For 8 = 0, 10, 80, 90 de- 
grees, the profile S(z) strikingly confirms IN phase coex- 
istence. We can clearly identify one region where S(z) 
is close to zero, corresponding to the isotropic phase, 
and another region where S(z) is closer to unity, corre- 
sponding to the nematic phase. Moreover, in the nematic 
phase, a(z) is roughly constant, and the plateau value 
closely follows the imposed anchoring angle 8. In other 



words, the constraint has the expected effect, namely to 
force the nematic phase to assume a specified tilt angle. 
Of course, in the isotropic phase, there is no preferred 
direction, and a(z) fluctuates randomly; one can show 
that the average should converge to 90(7r — 2)/7r « 32.7 
degrees. In contrast, for 8 = 45, 50 degrees, the scenario 
is completely different. Here, S(z) is roughly constant at 
S ~ 0.8, implying a single nematic phase along the entire 
z-direction. In addition, from the corresponding a(z), we 
see that the nematic is twisted: starting at z = 0, a(z) 
rotates smoothly from to 90 degrees, abruptly drop- 
ping back to again as one passes through the periodic 
boundary at z = 100. Clearly, this configuration does not 
reflect IN coexistence at all, but rather a twisted nematic 
phase with a surface defect. 

In light of Fig.[Tni it is clear that the free energy barrier 
for "in-between" tilt angles does not reflect the interfacial 
tension, and consequently Eq. (J2]) does not apply. For an- 
gles that are close to homogeneous and homeotropic an- 
choring, however, the IN scenario of Fig. [3] is confirmed. 
Therefore, for these angles, we may use Eq.([2]) to obtain 
the angle dependent interfacial tension j(8). The result is 
shown in Fig. [5J which reveals several trends. First of all, 
in contrast to v = 0, we now observe a profound variation 
of 7(0) with the imposed tilt angle. The interfacial ten- 
sion is smallest at 8 — 0, corresponding to homogeneous 
anchoring. We therefore expect unconstrained simula- 
tions, whereby 8 is not imposed but freely fluctuating, 
to mostly exhibit homogeneous anchoring. However, the 
data of Fig. [5] also suggest the presence of a shallow min- 
imum at 8 = 90, which corresponds to homeotropic an- 
choring. In other words, for v = 0.5, homeotropic anchor- 
ing appears to be metastable. Since the difference in in- 
terfacial tension between homogeneous and homeotropic 
anchoring is small, it is not a-priori clear which anchor- 
ing condition will actually prevail in an unconstrained 
simulation. 

We have therefore performed a number of uncon- 
strained simulations, i.e. without Eq.([3]), and measured 
the coexistence distribution In P(S). In addition, for each 
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simulation, we also recorded the n z component of the di- 
rector n as a function of S. In some cases, we found that 
the system selects 6 — 0, in which case n z drops to zero 
once nematic order sets in, but quite often also 9 = 90 is 
selected, in which case n z becomes close to unity. More 
precisely, using lattice dimensions L = 15 and D = 40, we 
performed 90 unconstrained simulations and found that 
metastable homeotropic anchoring (9 — 90) was selected 
27 times, i.e. in 30% of the cases. This finding is impor- 
tant because it shows that the order parameter distri- 
bution of the unconstrained simulation actually reflects 
a "weighted average" of both stable and metastable an- 
choring. This feature is illustrated in Fig.[9fa), which also 
includes In P(S) of the unconstrained simulation. As the 
figure shows, the free energy barrier of the unconstrained 
simulation is somewhere "in-between" homogeneous and 
homeotropic anchoring. 

Another important finding is that the unconstrained 
simulations reveal only stable anchoring (9 — 0), and 
metastable anchoring (9 — 90), while no other anchoring 
angles were observed. This result is consistent with "f(9) 
of Fig. [51 which indeed features just two minima. In other 
words, all "in-between" tilt angles are unstable. Systems 
in which the anchoring is held artificially fixed at such 
unstable angles, for example via the constraint of Eq.([3]), 
will experience an additional strain. For highly unsta- 
ble tilt angles, the strain is so strong, that it becomes 
favorable for the system to break-up the IN interfaces al- 
together, and form a twisted nematic. This is precisely 
the effect we observed for 9 = 45, 50 in Fig.[]I)J However, 
also for tilt angles close to the stable and metastable an- 
gle, we noticed that the strain manifests itself. In this 
case, a small shift in the transition inverse temperature 
e* can be detected. The effect is illustrated in Fig. \T[b), 
which shows e* as a function of the imposed tilt angle 9. 
For unstable tilt angles, e* is systematically larger com- 
pared to the stable and metastable angles. Of course, 
for the stable and metastable angles, which are the ex- 
perimentally relevant cases, one finds the same transition 
temperature again. Note also Fig.[7Ja), which shows that 
the effect for the spatially isotropic potential v = does 
not occur, as expected. 



E. Anchoring effects for v = —0.35 

For completeness, we also performed a number of sim- 
ulations using a negative value of v in Eq.Q, namely 
v = —0.35. Recall that for negative values, the side- 
side arrangement of Fig. 3] becomes energetically more 
favorable. Compared to v = 0.5, one might intuitively 
expect that this reverses the stable and metastable an- 
choring angles. The angle dependent interfacial tension 
indeed confirms this, see Fig. [5] We now observe that 
homeotropic anchoring yields the lowest interfacial ten- 
sion, i.e. is stable, while homogeneous anchoring appears 
to be metastable. In agreement with v = 0.5, we again 
measure a shift in e* when unstable anchoring angles are 
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FIG. 11: Coexistence distributions lnP(S') of Eq.(j4]|, with 
p = 10 and v — —0.35, using box dimensions L — 15 and 
D — 40. Shown are distributions for imposed tilt angles 9 = 
0, 90 degrees, as well as the unconstrained distribution that 
one obtains when the tilt angle is allowed to freely fluctuate. 



imposed, see Fig. [TJc). Interestingly, even though for 
v = —0.35 homeotropic anchoring yields the lowest in- 
terfacial tension, we observed that unconstrained simula- 
tions have difficulty "finding" this configuration. During 
a series of 95 unconstrained simulation runs, homeotropic 
anchoring was selected only 35 times, i.e. in 37% of 
the cases. In other words, the interfacial tension ex- 
tracted from lnP(S') in the unconstrained simulation, 
rather reflects the metastable anchoring condition, see 
Fig. HU The figure clearly shows that homeotropic an- 
choring (9 = 90) yields the lowest free energy barrier, 
while the barrier in the unconstrained distribution is sig- 
nificantly higher (and, in fact, rather closely resembles 
homogeneous anchoring) . From a computational point of 
view, the result of Fig. [TT] is important because it shows 
that simulations do not generally find the optimal an- 
choring angle by themselves, even in a relatively simple 
lattice model. 



IV. RESULTS: SOFT RODS 

Next, we investigate anchoring effects in an off-lattice 
fluid of soft rods. The rods are modeled as sphero- 
cylinders, of length I and width w. In this section, we 
set l/w = 10, and w will be the unit of length. The 
rods interact via a repulsive pair potential, whereby rod 
overlap is penalized with an energy cost of 2k^T. For 
more details about the model, the reader is referred to 
previous work [ID, The rods are simulated in the 

grand-canonical ensemble, i.e. at constant temperature 
T, chemical potential \x, and system volume V , while 
the number of rods in the system fluctuates. Again, we 
use an elongated simulation box V = L x L x D, with 
periodic boundary conditions. The simulations are per- 
formed using standard insertion/deletion moves [34j |. and 
the distribution lnP(S') is recorded, defined as the prob- 
ability to observe the nematic order parameter S, at the 
specified tilt angle 9. As before, 9 is imposed using the 
constraint of Eq. . For soft rods, we noticed that a sub- 
stantially larger threshold angle was needed to maintain 
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FIG. 12: Anchoring properties of soft rods at the IN transi- 
tion. Shown in (a) is the coexistence chemical potential fi* 
versus 0; in (b) the angle dependent interfacial tension "f(9) 
versus 0. Closed squares are raw simulation data; the curves 
serve to guide the eye. The horizontal line in (b) marks the 
interfacial tension of an unconstrained simulation, taken from 
previous work The simulations were performed using 

box dimensions L — 35 and D — 105. 



efficiency. Here, we used S = 2.5 degrees. Whereas in the 
thermotropic liquid crystal of Eq. dU) phase coexistence is 
achieved by tuning the inverse temperature e, here that 
role is played by the chemical potential fi. At the coex- 
istence chemical potential /i*, lnP(5) becomes bimodal: 
coexistence properties and interfacial tensions may then 
be extracted from the peak positions and heights, as in 
Fig.H 

The results of the soft rod simulations are summarized 
in Fig. [T2] Compared to the lattice simulations of Eq.((3]), 
the data reveal significant scatter. This indicates that 
soft rod simulations are demanding, and already close to 
the limit of what is currently tractable. Nevertheless, a 
number of trends emerge. According to Fig. [T27bh ^(9) 
increases monotonically with 9, with the minimum oc- 
curring at 9 = 0. In other words, soft rods favor homo- 
geneous anchoring, and the presence of metastable an- 
gles is unlikely. The data also show that the anchor- 
ing angle is a remarkably "soft" degree of freedom: the 
free energy cost of tilting the nematic director away from 
the IN interface is small. This is apparent from the co- 
existence chemical potential, see Fig. IT2T a). Note that 
Fig.[T27a) is the "analogue" of Fig.[7]for the lattice model 
of Eq.Q. For the lattice model, the coexistence inverse 
temperature increases profoundly away from the stable 
and metastable angles. This increase is a manifestation 
of the strain introduced into the system when unstable 
anchoring angles are imposed. In contrast, for soft rods, 
the coexistence chemical potential remains nearly con- 
stant over a wide range; only when 9 > 30 or so, does 
ix* begin to exhibit a pronounced 9 dependence. For soft 
rods, the anchoring angle can thus be varied around the 
stable direction over a fairly large range, without intro- 



ducing excessive strain into the system. This result is 
important for unconstrained simulations, where the an- 
choring angle is allowed to fluctuate freely. It is unlikely 
that such simulations would always reveal homogeneous 
anchoring. Rather, we expect a range of anchoring an- 
gles < 9 < 30 to be present. The horizontal line in 
Fig. fT2Tb) marks the interfacial tension obtained during 
an unconstrained simulation of soft rods (l9j , and indeed 
confirms this expectation. Even though the lowest in- 
terfacial tension is obtained at 9 — 0, the unconstrained 
simulation slightly exceeds this value. Instead, it rather 
reflects the average of 7(6*) over the range < 9 < 30 de- 
grees. Additional confirmation is obtained from simula- 
tion snapshots of unconstrained simulations, which reveal 
substantial fluctuations of the anchoring angle around 
the homogeneous direction. 



V. SUMMARY AND OUTLOOK 

In this paper, an alternative simulation approach to 
study anchoring effects at the IN interface in liquid crys- 
tals was described. In particular, we focused on the angle 
dependent interfacial tension "f(9), with 9 the anchoring 
or tilt angle. The proposed method is based on recent in- 
novations [H, [l!| where the order parameter distribution 
is used to extract interfacial properties. The new twist 
has been to introduce a constraint into the Hamiltonian, 
see Eq. ([3]) , which forces the nematic director to maintain 
a specified angle with respect to the a;y-plane. The idea 
is that, by using a simulation box that is elongated in 
the z-direction, IN interfaces will form in the xy-plane as 
well. The constraint then allows the anchoring angle 9 
to be fixed to some value of interest. 

At the same time, a new liquid crystal model was in- 
troduced. The model is defined on a lattice and exhibits 
a strong first-order IN transition. In addition, the pre- 
ferred anchoring (homogeneous, homeotropic, or neutral) 
can be tuned by means of a single parameter. Com- 
pared to more elaborate off-lattice models, such as rods 
or platelets, the lattice variant is considerably easier to 
simulate. In particular, equilibration is less problematic, 
and high-quality data are readily generated. Precisely 
this property was exploited to obtain j(9) for the lattice 
model, using the new method. Indeed, when anchoring 
effects are "switched-off " , by setting v = in Eq.©, 
7(0) becomes constant. In contrast, when v ^ 0, a pro- 
nounced 9 dependence is revealed. For these cases, only 
homogeneous and homeotropic anchoring were seen to be 
relevant. More precisely, for v > 0, homogeneous anchor- 
ing is stable, and homeotropic anchoring metastable. For 
v < 0, the trend is reversed. In other words, the preferred 
anchoring depends sensitively on the details of the inter- 
actions. Our results have also shown that, when unstable 
anchoring angles are imposed, the new method must be 
used with some care. In those cases, the simulations do 
not reveal IN coexistence, but rather a twisted nematic 
phase. Fortunately, when this happens, the method gives 
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a clear warning, in the form of a shift in the coexistence 
temperature. A somewhat surprising finding was that, 
even for the simple lattice model, simulations do not gen- 
erally find the "optimal" anchoring angle by themselves. 
Instead, when the nematic director is allowed to fluctuate 
freely, both stable and metastable anchoring are typically 
revealed. 

We have also applied the new method to obtain ^(9) for 
a fluid of soft rods. For soft rods, anchoring effects could 
also be identified, albeit that the data are significantly 
less accurate. The simulations reveal homogeneous an- 
choring to be stable, a finding which is consistent with 
most theoretical studies of hard rods. Interestingly, for 
soft rods, no metastable anchoring angle could be de- 
tected, which makes this model qualitatively very differ- 
ent from the lattice model of Eq.Q. It confirms, once 
again, that anchoring effects are extremely sensitive to 
the particle interactions. 



For the future, investigations of the capillary wave 
spectrum for the lattice model of Eq.(|3]) are planned. As 
mentioned in the Introduction, the spectrum is qualita- 
tively affected by the anchoring condition 0,@,0|- Since, 
in Eq.Q, the anchoring can be tuned using a single pa- 
rameter, and since the model is easy to simulate anyhow, 
such investigations should be worthwhile. A sound un- 
derstanding of the lattice model may well be a prereq- 
uisite before more complicated off-lattice simulations are 
attempted. 
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